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Abstract 


A  possible  countermunition  against  kinetic  energy  (KE)  rods  is  a  thin  square  metal  launched 
edp-on  so  that  an  edge  strikes  die  rod  first  followed  by  the  rest  of  the  plate.  A  properly  designed 
coil  gun  can  laimch  plates  edge-on  to  velocities  of  several  hundred  meters  per  second.  The 
proper  design  of  the  coil  gun  depends  on  knowing  the  mutual  inductance  between  the  coil  and 
die  plate,  which  was  calculated  by  assuming  that  the  current  distribution  in  plate  can  be  described 
by  a  polynomial  that  has  some  arbitrary  coefficients.  These  coefficients  were  determined  by 
equating  the  magnetic  induction  of  the  current  distribution  to  the  negative  of  the  applied 
magnetic  induction  of  die  launch  coil  at  selected  points  on  the  plate.  This  step  was  aided  by  the 
fact  that  all  the  integrals  arising  from  the  Biot-Savart  law  are  analytic.  In  the  next  step,  the 
mutual  inductance  was  calculated  from  the  current  distribution  and  the  applied  magnetic 
induction.  The  calculation  was  repeated  for  various  plate  positions  in  the  coil.  These  results 

were  used  to  design  the  power  supply  for  the  launch  coil  and  to  calculate  the  final  velocity  of  the 
plate. 
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1.  Introduction 


In  previous  work  in  the  area  of  plate-rod  interactions,  such  as  in  reactive  armor  [1]  and 
momentum  transfer  armor  (MTA)  [2, 3],  explosives  were  used  to  launch  the  plates  with  a  velocity 
in  the  direction  that  is  normal  to  the  plane  of  the  plate  (face-on  launch).  Launching  the  plate  with 
a  velocity  in  the  direction  that  is  in  the  plane  of  the  plate  (edge-on  launch)  was  not  considered, 
because  a  substantial  amount  of  explosive  materials  would  be  needed  behind  a  small  surface  area 
of  the  plate,  and  confin^ent  of  the  high-pressure  gases  would  be  difficult  There  is  some  evidence, 
however,  that  plates  launched  edge-on  may  be  more  effective  at  defeating  a  kinetic  energy  (KE)  rod. 
Both  Hull  code  calculations  [4, 5]  and  an  experiment  that  simulated  plates  flying  edge-on  by  firing 
a  yawed  projectile  against  a  stationary  target  plate  set  at  an  angle  [6]  suggest  a  preference  for  edge-on 
orientation.  This  was  confirmed  in  a  series  of  experiments  [7]  where  plates  were  launched, 
nonexplosively,  by  an  electromagnetic  (EM)  launcher  and  intercepted  a  subscale  rod.  There  are 
other  advantages  as  well.  A  plate  launched  edge-on  may  be  in  an  aerodynamically  preferable 
orientation  that  can  fly  the  distances  required  for  some  applications,  and  EM  launchers  do  not  have 
the  special  logistics  and  handling  difficulties  that  are  required  for  explosives.  Because  of  these 
advantages,  they  appear  to  be  a  viable  cation  and  warrant  continued  examination.  Accordingly,  this 
launch  option  for  KE  threats  is  being  studied  at  the  U.S.  Army  Research  Laboratory  (ARL). 

An  EM  launcher,  which  is  well  suited  for  edge-on  launch  of  plates,  is  the  reconnection  gun,  a 
type  of  coil  gun  invented  by  Cowan  [8].  In  this  launcher,  a  time  varying  magnetic  induction 
produced  in  an  external  laimch  coil  induces  a  current  in  the  plate  to  be  launched.  This  action  is 
similar  to  that  of  a  transformer  where  the  primary  winding,  the  external  launch  coil,  induces  a  current 
in  the  secondary  winding,  the  plate.  The  force  between  the  induced  current  and  the  external 
magnetic  induction  accelerates  the  plate  out  of  the  coil.  Cowan  and  colleagues  [9]  demonstrated  that 
high  velocities  can  be  achieved  by  using  a  number  of  external  coils  arranged  in-line  along  a  path. 
As  a  plate  was  passing  dirough  a  coil,  a  current  pulse  was  delivered  at  the  proper  time  to  accelerate 
the  plate  toward  the  next  coil.  This  resulted  in  the  launching  of  a  150-g  aluminum  plate  to  a  veloci^ 
of  1.0  km/s.  Although  this  multistage  laimcher  did  not  use  explosives,  its  use  in  nonenergetic 
momentum  transfer  armor  (NEMTA)  is  probably  not  practical  because  of  its  size  and  weight 
Therefore,  only  single-stage  launchers  are  considered  here. 
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2.  Single-Stage  CoU  Guns 


The  electrical  schematic  of  a  single-stage  coil  gun  [10]  (Rgure  1)  is  very  similar  to  a  series  LRC 
circuit  Unlike  the  classic  LRC  circuit,  tfie  inductance  associated  with  the  launcher  depends  on  the 
position  of  the  plate  within  the  coil.  This  dependence  is  due  to  the  distribution  of  the  magnetic 
induction  in  the  core  of  the  coil  and  around  the  plate.  The  resistor  in  this  figure  is  not  a  physical 
resistor.  It  represents  the  energy  losses  within  the  system  that  may  also  vary  with  time,  but  it  is 
usually  assumed  to  be  constant  for  the  purposes  of  modeling. 


SWITCH  RESISTOR 


Figure  1.  Schematic  of  a  Single-Stage  Coil  Gun. 


hi  designing  a  coil  gun  for  the  best  performance,  it  is  necessary  to  calculate  flie  current  through 
the  coil,  the  plate’s  velocity,  and  the  plate’s  position  in  flie  coil.  These  calculations  require  that  the 
inductance  of  the  coil,  and  its  gradient,  be  known  as  functions  of  the  plate’s  position.  The  coil’s 
inductance  is  necessary  to  calculate  the  coil’s  current,  and  the  inductance  gradient  is  necessary  to 
calculate  the  acceleration  of  the  plate.  Previously,  models  of  proposed  launch  coil  designs  were 
made  out  of  aluminum,  then  their  inductances  were  measured  with  a  Hewlett-Packani  model  4263 A 
LRC  meter.  This  meter  supplies  an  alternating  voltage  to  the  coil  and  analyzes  the  amplitude  and 
phase  of  the  current  tiirough  tiie  coil  to  determine  its  inductance.  The  results  are  digitally  diq)layed 
on  the  meter  and  recorded  for  each  plate  position.  These  data  are  then  used  to  estimate  the 
inductance  gradients.  This  study  was  tedious  because  of  the  number  of  coils  that  had  to  be 
constmcted  and  measured,  ff  it  were  possible  to  easily  calculate  the  inductance  and  the  gradient  of 
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these  coils,  then  a  more  complete  study  of  various  coils  with  different  design  parameters  and 
geometries  could  have  been  performed.  This  report  presents  a  method  to  calculate  these  quantities. 

To  illustrate  the  change  in  die  distribution  of  the  magnetic  induction  when  a  conductor  is  present. 
Figure  2  shows  the  result  of  a  simple  two-dimensional  calculation  [10].  The  shaded  area  represents 
an  aluminum  plate  in  an  oscillating  magnetic  induction,  produced  by  two  in&iite  current  sheets, 
located  at  the  thick  solid  lines  to  the  left  and  to  the  right  of  the  aluminum  plate.  Current  in  these 
sheets  is  equal  in  magnitude,  but  opposite  in  direction  at  all  times.  The  total  magnetic  induction 
produced  by  these  current  sheets  and  the  eddy  current  in  the  plate  is  shown  as  streamlines,  located 
between  the  current  sheets.  The  direction  of  the  magnetic  induction  is  tangent  to  diese  lines,  and  the 
magnitude  is  inversely  proportional  to  the  distance  between  neighboring  lines.  Because  the  magnetic 
induction  produced  by  die  eddy  current  tends  to  be  equal,  but  opposite  to  the  magnetic  induction  of 
the  current  sheets,  their  sum  is  small  inside  the  aluminum  plate,  except  at  the  left  and  right  edges 
where  the  streamlines  may  ent»  the  plate.  Thus,  the  streamlines  are  bent  around  the  aluminum  plate. 
This  bending  of  the  streamlines  concentrates  diem  in  the  gap  between  the  back  edge  of  the  plate  and 
the  left  current  sheet,  indicating  an  area  where  there  is  a  large  magnetic  induction  and  force  on  the 
plate.  This  force  pushes  the  plate  through  the  slot  in  the  ciurent  sheets  to  the  right  and  out  of  the 
coil.  As  the  plate  moves  out,  the  streamlines  straighten  out  and  eventually  become  evenly  spaced 
vertical  lines,  when  the  plate  is  completely  out  of  the  coil.  The  shape  of  die  streamlines  and  the 
force  on  the  plate  can  be  easily  calculated  for  this  two-dimensional  geometry.  Unfortunately,  two- 
dimensional  calculations  are  not  adequate,  as  has  been  demonstrated  by  another  more  detailed  two- 
dimensional  calculation  [11],  where  the  force  on  the  plate  for  a  transient  current  in  the  sheets  and 
the  heating  effects  in  the  plate  were  included.  The  plate  velocity  was  calculated  to  be  140  m/s,  but 
the  experimental  velocity  was  only  92  m/s.  There  was  better  agreement  when  an  “effective  3-D 
correction  factor”  was  introduced. 

3.  Eddy  Currents 

Three-dimensional  eddy  current  calculations  [12]  are  important  in  electrical  engineering,  and 
many  publications  are  available  on  diis  subject  Most  of  these  calculations  can  be  divided  into  two 
broad  categories.  Hie  first  category  is  the  quasi-stationary  [13]  problem,  where  the  magnetic  field 
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Figure  2.  Magnetic  Streamlines  Around  an  Aluminum  Plate. 


is  harmonically  oscillating  at  a  constant  ftequency  and  amplitude.  The  second  category  is  the 
transient  [14]  problem,  where  die  magnetic  field  pulse  begins  at  zero,  varies  with  time,  and  returns 
to  zero  at  some  later  time.  Some  publications  present  computer  programs  that  can  solve  either 
category  of  problems  for  many  geometries.  These  noncommercial  computer  programs,  however,  are 
complicated  and  difficult  to  use.  Furthermore,  the  utility  of  these  programs  is  limited,  since  it  is  not 
always  possible  to  make  simplifying  assumptions  pertinent  to  the  problem  at  hand.  As  an  example, 
a  general  treatment  of  the  subject  would  include  the  conductivity  of  the  plate  and  the  time  variation 
of  the  applied  magnetic  field,  but  these  features  may  not  be  necessary  for  this  application. 

These  features  were  not  included  in  a  preliminary  calculation  [15],  where  it  was  assumed  that 
the  eddy  current  of  the  plate  was  distributed  only  on  the  surface  of  the  plate  and  the  magnetic  field 
was  static.  To  model  the  eddy  current,  the  plate  was  covered  on  the  top  and  the  bottom  by  an  array 
of  filamaitary  rectangular  loops.  To  model  the  condition  that  the  magnetic  induction  produced  by 


the  eddy  current  is  equal  and  opposite  to  the  coil's  magnetic  induction  in  the  plate,  the  current  in  each 
loop  was  chosen  to  make  the  total  magnetic  induction  zero  at  the  center  of  each  loop.  Once  the 
rnm»nt  in  each  loop  was  calculated,  the  force  on  each  loop  was  found  and  summed  to  give  the  total 
force  on  die  plate  that  is  dependent  on  the  inductance  gradient  [10], 


Fix)  = 


L'jx) 


(1) 


where  Fix)  is  the  force  on  the  plate,  I  is  the  current  in  the  coil,  and  L'ix)  is  the  inductance  gradient 
of  the  coil.  This  calculation  was  repeated  for  a  number  of  plate  positions,  and  the  results  were 
compared  with  the  inductance  gradient  as  determined  from  the  inductance-meter  measurements. 
This  procedure  was  repeated  for  a  number  of  plate  dimensions,  plate  materials,  and  coil  designs.  The 
calculated  inductance  gradients  compared  favorably  with  the  measured  inductance  gradients  in  all 
cases,  provided  that  the  plate  was  a  good  conductor  and  had  a  thickness  larger  than  the  skin  depth. 


These  results  dmonstrate  that  a  current  distribution  on  the  surface  may  be  a  good  approximation 
for  the  actual  current  distribution  within  the  plate.  Indeed,  the  actual  eddy  current  distributions 
found  in  other  simpler  problems  are  concentrated  within  some  depdi  from  the  surface,  ff  this  depdi 
is  small,  compared  to  the  dimensions  of  the  plate,  dien  this  approximation  can  be  made.  The  exact 
definition  for  the  current  depth  is  not  known  for  a  rectangular  plate,  but  it  can  be  estimated  by 
considering  a  simpler  problem.  Let  a  conductor  occupy  all  the  positive  x  space  in  the  three- 
dimensional  Cartesian  coordinate  firame  and  the  rest  of  the  space  be  occupied  by  a  harmonically 
oscillating  magnetic  induction  in  the  y  direction,  which  induces  an  eddy  current  in  the  conductor. 
The  analytical  solution  [16]  to  fliis  problem  diows  tiiat  the  amplitudes  of  tiie  magnetic  induction  and 
the  eddy  current  decay  exponentially  wifli  distance  from  die  surface,  with  a  decay  lengdi  or  “the  skin 
depth” 


6  = 


poo) 


(2) 


where  p  is  the  permeability  of  the  plate,  w  is  the  angular  frequency  of  the  magnetic  induction,  and 
o  is  the  conductivity  of  the  plate.  In  the  present  launching  systems,  die  fi^quency  is  typically  1  kHz, 
(0  =  6282  rad/s,  and  the  plates  are  aluminum,  o  =  3.54  x  10’  S/m.  Consequaitly,  the  skin  depth  is 
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2.7  mm.  This  dq)th  is  small,  compared  to  the  typical  length  and  width  of  the  plate  which  are  greater 
than  0.10  m,  but  it  may  not  be  small  when  compared  to  the  typical  plate  thickness  of  6  mm.  The 
orientation  of  the  plate  relative  to  the  magnetic  induction  (Figiu®  2)  illustrates  that  the  dimensions 
to  be  compared  to  the  skin  depth  are  die  length  or  width,  and  not  the  diickness.  the  skin  depth 
were  zero,  in  this  figure,  the  streamlines  that  enter  the  plate  at  the  edges  would  be  completely 
expelled.  Other  than  this,  the  shape  of  the  rest  of  the  streamlines,  elsewhere,  would  change  very 
little. 

A  zero  skin  depth  will  result  if  either  the  angular  fi»quency  6)  approaches  infinity,  “the  higb- 
frequency  limit,”  or  if  the  conductivity  of  the  material  o  approaches  infinity,  “the  super-conductor 
limit.”  In  the  super-conducting  limit,  the  angular  frequency  can  be  zero,  which  means  that  the 
magnetic  induction  field  can  be  static.  This  also  makes  the  problem  easier  to  solve  but  raises  an 
interesting  question;  How  can  eddy  currents  be  produced  by  a  static  magnetic  induction  field?  Eddy 
currents  are  induced  in  a  super  conductor  in  a  different  manner.  Consider  a  super  conductor  with 
no  eddy  currents  located  in  a  field-free  region,  far  fi*om  the  core.  Now,  bring  the  super  conductor 
into  die  core  where  diere  is  a  magnetic  induction.  This  motion  produces  the  time  varying  magnetic 
induction  that  induces  the  eddy  current  Thus,  the  time  variance  of  the  magnetic  induction  is  implied 
by  the  condition  that  the  magnetic  induction,  inside  the  super  conductor,  is  zero,  which  means  that 
it  had  to  be  in  a  field-free  region  at  some  time.  Furthermore,  because  work  must  be  performed  to 
establish  the  eddy  currents  themselves,  die  super  conductor  must  be  pushed  into  the  core  and  it  will 
be  expelled  out  of  the  core  if  it  were  released.  Thus,  a  static  magnetic  field  can  launch  a  super¬ 
conducting  plate,  which  is  the  “Meissner  effect.” 

4.  Eddy  Current  Calculations 

The  current  loops  in  the  preliminary  calculation  are  now  replaced  by  a  continuous  current 
distribution,  and  it  is  assumed  that  the  rectangular  plate  has  no  thickness.  The  magnetic  induction 
produced  by  an  eddy  current  in  this  thin  plate  is  given  by  the  Biot-Savart  law: 
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(3) 


47cp^  J 


fd.'fdy 

•fl  -b 


,  j (x'y)x(x-r) 
lie-  jc'P 


where  a  is  half  the  length  and  b  is  half  the  height  of  the  plate.  If  the  plate  has  infinite  conductivity, 
thai  this  magnetic  induction,  ( Jc),  must  be  equal  and  opposite  to  the  coil’s  magnetic  induction,  (x), 

everywhere  on  the  plate  so  that  the  total  is  zero  everywhere  in  this  area.  This  cannot  be  done  if  tiie 
coil’s  magnetic  induction  has  a  component  in  the  x  or  y  direction  on  the  plate.  When  x  -  x '  is  on 
thex  -  y  planft,  the  emss  product  in  the  numerator  has  just  one  component  in  the  z  direction,  and  so  B^ 
must  also  be  in  this  direction.  Thus,  there  cannot  be  an  x  or  a  y  component  to  cancel  out  flie  like 
components  of  the  coil’s  magnetic  induction.  Fortunately,  the  plates  are  positioned  so  that  these 
components  of  the  coil’s  magnetic  induction  are  small  compared  to  the  z  component  This 
complication  would  not  occur  if  the  plate  had  a  thickness,  and  the  current  density  on  all  six  faces  of 
the  plate  were  included. 

The  problem  to  be  solved  here  is  more  complicated  than  the  usual  magnetostatic  problem  where 
the  current  density  is  given  first  and  one  must  then  find  the  magnetic  induction  everywhere.  In  this 
case,  the  coil’s  or  the  eddy  current’s  magnetic  induction  is  given  first  and  one  must  then  find  the 
current  density.  One  way  to  solve  this  problem  is  to  assume  that  the  current  density  may  be 
approximated  as  a  sum  of  terms  that  are  products  of  a  coefficient  and  a  member  of  a  set  of  linearly 
independent  functions.  As  an  example,  the  x  component  of  the  current  distribution  is  assumed  to 
be 

C4) 

m,ii=0 

and  the  y  component  of  the  current  distribution  is  assumed  to  be 

Jy(.x\y')  =  E  (5) 

/,*=o 
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The  set  of  linearly  independent  functions  in  this  case  consists  of  products  of  x,  raised  to  an 
integer  power,  and  y,  raised  to  an  integer  powCT.  When  these  functions  are  substituted  into 
equation  (3),  the  resulting  integrals  are  analytic,  and  calculations  are  simplified  (see  the  Appendix). 

These  coefficients  must  satisfy  several  conditions.  First,  the  current  density  that  they  describe 
must  conserve  charge  everywhere  on  the  plate,  or 


a/,(x',y')  dJ(x\y) 
dx'  ^  dy' 


=  0, 


(6) 


for  any  point  on  the  plate.  Second,  the  ciurent  density  must  not  cross  the  edges,  which  means  that 
no  cunent  may  enter  or  leave  the  isolated  plate.  Thus,  (x'  =  a,  y')  must  be  zero  for  all  points  on 
the  right  edge  and  (pc'  =  -a,  y')  must  be  zero  for  all  points  on  the  left  edge.  Similar  conditions 
apply  to  the  top  and  bottom  edges,  (x',  y'  =  ±b)  =  0.  The  current  density  vector,  however,  can  be 
parallel  to  an  edge.  These  conditions  make  some  of  die  coefficients  dependent  on  each  other; 
choosing  a  value  for  one  coefficient  will  change  die  value  of  another.  It  is  possible  to  derive  the 
relations  between  the  coefficients,  but  there  is  another  easier  method  to  generate  a  physically 
meaningful  current  distribution.  In  this  method,  die  current  distribution  is  generated  ftrom  another 
linear  combination  of  the  basic  functions: 


u 


Tix',y')  =  (a^-x'^)  (b^-y'^)  T.^x'*  y'\ 


i.J‘0 


where 


^  ax 


and 


JAx',y') 


dT(x',y') 

dy' 


(7) 


(8) 


(9) 
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I  and  J  in  equation  (7)  are  the  maxiniuni  powers  of  x'  and  y'  for  the  sununation.  Their  exact  values 
depend  on  the  desired  degree  of  approximation.  For  example,  I  and  J  were  equal  to  9  in  the  results 
that  will  be  presented  later. 

To  illustrate  how  equation  (7)  meets  the  conditions  for  a  physical  current  distribution,  consider 
a  current  distribution  described  by  a  single  term, 


T.  .  (x',y)  =  (a^  -  x'^)  (b^  -  y'2)  y'^  T.  j, 


(10) 


where 


dT  (x' 


(11) 


and 


Jy  (Jc',/)  =  =  -T  (1,2- y'2'^  x'*'^y^iia^-(i  +  2)x'^). 

y  *•' 


(12) 


This  current  distribution  conserves  charge  and  satisfies  the  boundary  conditions  at  the  edges  for  all 
values  for  i  and  j.  Now  that  each  term  in  equation  (7)  generates  a  physical  current  density 
distribution  and  their  sum  will  also  generate  a  physical  current  density  distribution,  the  coefficients 
Tij  are  independent  from  each  odier.  When  equations  (1 1)  and  (12)  are  substituted  into  the  Biot- 
Savart  law,  equation  (3),  the  magnetic  induction  everywhere  for  this  current  distribution  is 


B^(x,y,z\j 


/  (y - yO  (g ^ - x'^) x'^y'^' ^(Jb^-U+2) y'^) 
4n  ((x-x'f  +  iy-y'y+z^)^^^ 


T  U  ", 
4it 


jdx'fdy 

-a  -b 


,  (x-x')(b^-y^)x'^-^y^iia^-ii  +  2)x'^) 


ax-x'^Hy-yr-^z^) 


>\2^,2\3I2 


(13) 


As  stated  before,  all  these  integrals  are  analytic  and  procedures  for  their  evaluation  are  given  in 
the  Appendix.  The  integrals  in  equation  (13)  are  collected  and  defined  as 
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-a  -b 


Ay-y')(a^-x'^)x'‘yj-^Ub^-a-^2)y'^) 

((x-x')^+(y-y')^+Z^)^'^ 


+ 


,ix-x')(b^-y'^)x'*-^y\ia^-ii+2)x''^) 
((x-x')^  +  (y-y')^+Z 


(14) 


Thus,  if  all  the  coefficients  were  known,  then  the  z  component  of  the  magnetic  induction  of  Ae  eddy 
current  everywhere  would  be 


u 

Uj-0 


B,(x.y,z)  =  £  T,j  I,(x.y.z),,. 


(15) 


Similar  expressions  may  be  written  for  the  x  and  y  components,  but  they  are  not  needed  because  they 
are  zero  when  z  =  0.  By  choosing  a  set  of  field  points  on  the  plate,  it  is  possible  to  find  values  for 
the  coefficients  so  that  the  magnetic  induction  produced  by  the  eddy  current,  equation  (15),  is  equal 
and  opposite  to  the  coil's  magnetic  induction  at  these  points.  As  an  example,  assuming  that  /=  1  and 
/  =  1,  there  are  four  unknown  coefficients,  Too,  Tq,!,  Tj  o,  and  r,,i,  that  are  determined  by  four 
simultaneous  equations.  These  equations  result  from  making  (x,y,0)  =  (x,y,0)  at  four  points: 


5z(x„y„0)  =  -B,(x„y„0)  =  ro.o4(x„y„0)o^  +  To,,  /,(x,.y„0)o.,  +  /,(x„y„0),^  +  T,.,  /,(x,.y„0),.„ 

5,(x2,y2,0)  =  -B,(x2,y2,0)  =  To.o4(x2,y2,0)o.o+  To.,  4(x2,y2,0)o,,  +  r,.o  4(x2,y2,0)i,o  +  Ti.,  /,(x2,y2,0)i.„ 

5*(x3,y3,0)  =  -B,(x3,y3,0)  =  T’o.o  /x(x3.y3.0)o,o  +  To.,  4(x3,y3,0)o.,  +  /,(x3,y3,0),^  +  Ji.,  4(x3,y3,0)i.„ 

^z(X4.y4.0)  =  -Bc(X4.y4.0)  =  ro.o  /x(X4,y4.0)o,o  +  To.,  /x(X4,y4,0)o,,  +  /x(X4,y4,0),^  +  r,.,  /x(X4,y4,0)i.,.  (16) 


The  exact  positions  of  the  points  are  arbitrary,  as  long  as  they  are  not  on  the  plate's  edge.  Choosing 
the  points  on  a  uniform  grid,  such  as  the  one  shown  in  Figure  3,  an  example  for  a  higher  order  of 
approximation  (1  =  7  andJ  =  7),  works  well.  There  are  eight  rows  and  columns,  because  the  indices 
i  andj  start  at  0  and  continue  up  to  and  including  7,  in  this  illustration.  The  resulting  set  of  equations 
can  be  solved  by  using  standard  methods  in  linear  algebra  to  find  the  coefficients.  Because  these 
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Figure  3.  Grid  Points  Where  the  Total  Magnetic  Induction  Is  Zero. 

points  are  fixed  on  die  plate,  the  integrals  do  not  change  when  the  plate  is  moved  away  from  the  coil. 
Thus,  the  matrix,  whose  elements  are  the  integrals,  needs  to  be  evaluated  and  inverted  only  once. 
The  componaits  of  the  vector  that  must  be  multiplied  by  the  inverse  matrix,  however,  are  negative 
of  the  coil’s  magnetic  induction  at  each  point  Because  these  components  do  depend  on  the  relative 
position  between  the  coil  and  the  plate,  tiiey  must  be  reevaluated  when  the  plate’s  position  is 
changed. 

The  launch  coil  that  will  be  used  in  die  following  calculations  was  constructed  fix)m  two  copper- 
beryllium  alloy  plates  15.24  x  15.24  x  0.64  cm.  Each  plate  was  milled  into  a  square  helix  by  a 
0.64-cm-diameter  end  mill.  The  milling  pattern  left  a  conductor  with  cross-sectional  dimensions  of 
0.64  X  0.64  cm  in  a  square  helix,  with  five  complete  turns.  The  square  helixes  were  mounted  parallel 
to  each  other,  separated  by  5.08  cm.  The  magnetic  induction  of  the  coil  at  the  grid  points  was 
estimated  by  replacing  the  conductors  with  a  current-carrying  filament,  located  at  the  geometric 
center  of  their  cross  sections.  Although  current  is  actually  distributed  over  the  cross  section  of  the 
bars,  the  magnetic  induction  of  a  distributed  current  approaches  that  of  a  current  filament  at  some 
distance  from  the  bar.  Thus,  the  magnetic  induction  of  the  current-carrying  filaments  should  be 
approximately  the  same  as  the  coil,  provided  that  the  plate  is  not  too  close  to  a  coil’s  conductor. 
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In  Figure  4,  the  center  of  a  15-  x  15-cm  plate  was  positioned  at  the  center  of  the  coil.  The  lines 
within  the  plate  are  the  streamlines  for  a  current  density  that  cancels  the  coil's  magnetic  field  at 
100  grid  points,  7  =  9  and  J  =  9.  These  streamlines  have  the  same  properties  as  the  magnetic 
streamlines  in  Figure  2.  The  direction  of  the  current  density  is  tangent  to  the  line,  and  the  magnitude 
is  inversely  proportional  to  the  distance  between  neighboring  lines.  The  same  spacing  of  the 
streamlines  at  each  edge  of  the  plate  indicates  that  the  magnitudes  of  the  current  density  are  about 
the  same  at  each  edge.  Because  the  magnitudes  of  the  magnetic  induction  are  also  about  equal  at  the 
edges,  the  magnitudes  of  the  force  on  each  edge  are  about  equal.  These  forces,  however,  are  directed 
toward  the  center  of  die  plate,  and  the  plate  will  not  be  accelerated  out  of  the  coil  when  placed  in  this 
position.  In  contrast.  Figure  5  shows  the  streamlines  when  the  plate  is  half  way  out  of  the  coil.  The 
streamlines  are  now  concentrated  along  the  back  edge  of  the  plate,  indicating  an  area  where  there  is 
a  large  eddy  current  and  force  on  the  plate,  which  accelerates  the  plate  out  of  the  coil. 


Figure  4.  Current  Streamlines  for  a  Plate  at  the  Center  of  the  Coil. 
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Figure  5.  Current  Streamlines  for  a  Plate  Half  Way  Out  of  the  Coil. 

Eddy  currents  in  a  plate  were  evident  when  the  surface  temperature  distribution  of  a  plate  was 
observed  just  after  it  was  launched  in  an  earlier  experiment  [10].  The  launch  coil  was  made  from 
a  10-  X  10-cm  square  tube  of  aluminum  23-cm  long  with  3-mm-thick  walls.  Its  sides  were  slotted 
so  that  the  remaining  aluminum  formed  a  square  helical  coil,  with  nine  turns.  The  coil  was 
supported  on  the  inside  and  outside  by  a  stack  of  5-cm-thick  GIO  fiberglass  rectangles.  The  pacing 
between  one  pair  of  windings  was  widaied  in  order  to  permit  insertion  of  die  metal  plate  into  its  core 
from  the  side.  An  aluminum  plate  was  fabricated  with  a  nose  section  that  held  a  nail  oriented  along 

the  direction  of  plate  motion.  Ih  a  low-velocity  launch,  the  nail  was  driven  into  a  plywood  barrier 

* 

by  die  plate.  In  this  manner,  the  assembly  was  captured  for  viewing  by  an  infrared  (IR)  video  camera 
(8-12  pm)  immediately  after  latmch.  The  aluminum  plate  had  a  thick  anodized  coating  to  increase 
the  emissivity  of  the  surface  for  IR  radiation.  The  IR  video  image  in  Figure  6a  shows  that  the  top 
surface  of  the  plate  was  heated  along  its  trailing  and  side  edges.  By  obsaving  the  captured  plate  for 
some  time  with  the  IR  video  camera,  we  found  that  the  time  for  the  heat  to  diffuse  from  the  edges 
was  long,  compared  to  the  time  to  launch  the  plate.  Thus,  it  was  concluded  that  there  was  negligible 
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Figure  6a.  IR  Video  Image. 


heat  diffusion  during  launch.  Assuming  no  heat  diffusion,  and  asfaiTning  that  the  conductivity  of  the 
material  is  independent  of  the  temperature,  Ihe  heating  rate  is  then  proportional  to  the  square  of  die 
magnitude  of  the  current  density.  The  current  density  was  calculated  for  the  plate  at  its  initial 
position  inside  the  coil,  because  the  plate  remained  close  to  the  initial  position  during  the  entire 
current  pulse  for  this  low-velocity  shot.  The  shape  of  die  contour  lines  (Figure  6b),  for  the  heating 
rates  and  die  IR  video  image,  shows  that  the  side  edges  are  heated  by  large  eddy  currents  generated 
by  the  proximity  of  the  coil  windings.  No  attempt  was  made  to  correlate  die  heating  rate  with  the 
observed  intensity  of  the  IR  radiation. 

5.  Plate-Coil  Mutual  Inductance 

After  the  current  density  is  calculated,  the  mutual  inductance  can  be  found  from  the  general 
expression  for  the  energy  to  establish  a  magnetic  field. 
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Figure  6b.  Heating  Rate  Contour  Lines. 

=  O’) 


where  7g  is  the  current  density  of  the  coil,  is  the  plate's  current  density  or  die  eddy  current,  and  A 
is  the  total  magnetic  vector  potential.  The  first  integral  is  taken  over  the  volume  of  the  coil,  dv,.,  and 
the  second  is  taken  over  die  volume  of  die  plate,  dv^  When  the  plate  is  a  super  conductor,  the  total 
magnetic  induction  is  zero  everywhere  within  the  volume  of  the  plate,  and  the  magnetic  vector 
potratial  can  be  expressed  as  a  gradient  of  a  scalar  function,  A  =  grad  x.  Substituting  this  into  the 
last  integral,  it  is  possible  to  prove  that  the  integral  is  zero.  Thus,  the  energy  of  the  magnetic  field 
is  just  the  first  integral,  provided  that  the  plate  is  a  super  conductor.  The  next  step  is  to  substitute 
the  sum  for  the  total  magnetic  vector  potential. 


if 


dv^J, 


(18) 
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IS  the  vector  potential  produced  by  the  current  density  in  the  coil,  and  A^  is  the  vector  potential 
produced  by  the  eddy  current  in  the  plate.  It  is  now  possible  to  replace  each  term  in  this  equation 
with  an  equivalent  inductive  element.  First,  start  with  the  total  magnetic  energy  W.  If  the  plate  is 
placed  at  some  position  in  the  coil  Xq,  the  coil  has  some  inductance  L  (Xg),  and  if  a  current  /  were 
flowing  through  the  coil,  the  total  energy  stored  in  the  coil  would  be 


W  = 


2 


(19) 


The  first  integral  on  the  right,  in  equation  (18),  represents  the  magnetic  energy  stored  in  the  coil  itself 
because  it  depends  only  on  the  coil’s  current  density  and  geometry.  Thus,  this  integral  can  be 
replaced  with  a  similar  expression. 


(20) 


whereL,  is  the  self  inductance  of  the  coil.  Substituting  this  into  equation  (18),  the  inductance  of  the 
coil  becomes 


=  L,  + 


(21) 


The  integral  in  this  equation  is  not  in  the  most  convenient  form,  because  the  vector  potential  of  the 
eddy  current  is  used.  Although  it  can  be  calculated,  it  is  possible  to  rewrite  the  integral  in  terms  of 
quantities  that  have  already  been  calculated  in  the  course  of  finding  the  eddy  current.  By  using  the 
definition  of  the  magnetic  vector  potential  in  tarms  of  current  density  and  some  vector  identities,  the 
following  identities  can  be  proven  for  a  coil  and  plate  with  any  solid  shape,  each  having  a  current 
density  distribution: 


(22) 
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where  J  =  curl  T.  The  vector  T  must  be  either  zero  or  normal  to  the  surface  for  all  points  on  the 
surface.  These  general  identities  and  conditions  are  slightly  modified  for  the  present  application. 
The  vector  7  has  just  one  component  in  the  z  direction,  whose  curl  reproduces  equations  (8)  and  (9), 
and  it  is  zero  at  the  plate’s  edges.  In  addition,  only  the  z  component  of  die  coil’s  magnetic  induction 
is  used.  Thus,  the  inductance  of  the  coil  in  this  application  is 

L  \jdxjdyB^  {x,y)  T(x,y).  (23) 

^  -a  -b 


Because  the  coil’s  magnetic  induction  and  the  eddy  currents  are  proportional  to  the  coil’s  current  /, 
the  inductance  of  the  coil  will  be  independent  of  /.  Thus,  L  (x^)  will  depend  only  on  the  geometry. 

A  procedure  to  estimate  the  self-inductance  of  the  coil  L,  has  been  developed  and  presented  in 
detail  elsewhere  [15],  but  it  will  be  quickly  described  here.  This  [Hxx^ure  first  divides  the  winding 
of  the  coil  into  a  number  of  rectangular  bars  with  different  lengths  and  orientations.  The  self¬ 
inductance  of  the  coil  is  then  a  sum  of  the  self-inductance  of  each  bar,  and  the  mutual  inductance 
between  each  pair  of  bars.  The  self-inductance  of  a  bar  is  approximated  by  the  geometric  mean 
distance  formula  [17].  The  mutual  inductance  between  a  pair  of  bars  is  approximated  by  replacing 
each  bar  with  a  filament,  located  at  the  center  of  die  cross  section  of  the  bar,  and  using  the  exact 
formula  for  the  mutual  inductance  between  two  filaments  [17].  Using  these  approximations,  and 
taking  the  appropriate  sum,  the  self-inductance  of  the  launch  coil  being  considered  here  was 
calculated  to  be  5.90  pH.  Using  5.90  pH  for  the  self-inductance  of  the  coil  and  calculating  the 
integral  in  equation  (23)  for  various  plate  positions,  the  total  inductance  of  the  coil  as  a  function,  die 
plate’s  position  is  shown  as  the  upper  curve  in  Figure  7.  The  lower  curve  is  the  inductance,  when 
the  measured  self-inductance  of  the  coil  5.49  pH  is  used  instead.  The  diamonds  in  this  figure  are 
the  measured  inductance  of  the  coil  for  various  plate  positions. 
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6.  Discussion 


The  agreement  between  the  lower  curve  and  die  experimental  measurements  in  Figure  7  shows 
that  this  procedure  is  capable  of  estimating  the  mutual  inductance  between  the  coil  and  the  plate. 
It  also  shows  that  the  procedure  for  calculating  the  self-inductance  of  the  coil  gives  a  reasonable 
approximation.  Therefore,  it  is  now  possible  to  calculate  the  inductance  of  a  coil  as  a  function  of 
the  position  of  a  rectangular  plate  widiin  it  This  permits  study  of  various  coils  to  identify  what 
design  features  will  improve  the  efficiency  of  the  single-stage  coil  gun.  The  results  of  diis  study  will 
be  presented  in  future  publications. 


Figure  7.  Calculated  Coil  Inductance  (Solid  Lines)  and  Measured  Inductance  (Diamonds). 

Ii  practice,  it  is  found  that  the  series  expansion  of  the  cunent  streamline  function,  equation  (7), 
is  semiconveigent  When  the  mutual  inductance  of  the  plate  and  the  coil  for  a  given  plate  position 
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is  calculated  for  increasing  I  and  J  in  equation  (7),  the  mutual  inductance  seems  to  converge  to  a 
limit,  but  diverges  when  they  are  made  too  large.  This  divergence  may  be  caused  by  an  instability 
in  the  analysis.  Therefore,  I  and  J  should  be  kept  small  enough  to  avoid  the  instabilities  but  large 
enough  to  get  reasonable  results.  By  repeating  the  calculation  for  different  ‘T’s  and  ‘T’s  and 
observing  when  the  mutual  inductances  start  to  diverge,  it  is  possible  to  decide  what  their  values 
should  be.  This  instability  may  be  indicated  by  the  observations  that  the  coefficients  for  the  higher 
order  terms  do  not  approach  a  limit  with  increasing  I  and  /,  and  the  absolute  values  of  the 
coefficients  increase  as  the  order  of  their  terms  increases.  These  features  are  contrary  to  a 
converging  series  and  make  it  difficult  or  impossible  to  assign  an  error  to  the  results. 

The  reason  for  this  semiconvergence  may  be  illustrated  by  applying  this  method  to  a  simple 
analytic  problem.  Let  the  plate  in  Figure  3  be  an  infinitely  long  ribbon  along  the  x-axis,  with  a  half 
width  b,  and  assume  that  there  is  a  uniform  magnetic  field  in  the  z  direction  at  large  distances  away 
from  the  ribbon.  It  can  be  shown  analytically  that  die  current  density  in  the  ribbon  is 


(24) 


where  Bg  is  the  magnitude  of  the  magnetic  induction.  Because  the  current  density  is  infinite  at  the 
edges  of  the  ribbon,  it  is  quite  possible  that  there  are  infinite  current  densities  at  some  of  the  edges 
of  the  rectangular  plate.  If  this  is  so,  polynomials  of  finite  order  are  a  poor  approximation  for  the 
current  density  near  these  edges,  since  they  are  finite  everywhere  on  the  plate.  Because  the  distance 
between  these  edges  and  the  nearest  grid  points  deo^ases  as  7  and  J  are  increased,  it  is  possible  that 
the  higher  order  approximations  for  the  current  density  become  unstable.  Had  the  plate,  however, 
had  a  finite  thickness  and  conductivity,  thai  the  current  density  should  be  finite  everywhere.  These 
polynomials  may  then  be  a  very  good  approximation  for  the  current  density.  Thus,  this  formalism 
is  now  being  extended  for  this  case,  and  will  be  the  subject  of  future  reports. 
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7.  Conclusion 


This  calculation  of  the  eddy  currents  and  the  mutual  inductance  depends  on  three  assumptions: 
(1)  the  skin  depth  is  zero,  (2)  the  plate  has  zero  thickness,  and  (3)  flie  eddy-current  density 
distribution  in  the  plate  can  be  described  by  a  polynomial.  The  first  assumption  was  made  in  a 
preliminary  calculation  [15]  that  produced  very  good  agreement  with  experimental  measurements. 
These  results  indicated  that  the  skin  depth  may  be  small,  compared  to  some  dimension  of  the  plate. 
After  estimating  the  skin  depth  and  considering  the  distribution  of  the  magnetic  induction  around 
the  plate,  it  was  concluded  that  the  skin  depth  should  be  compared  with  the  lengtii  or  the  width  of 
the  plate,  and  not  with  its  thickness.  This  calculation  also  demonstrated  that  there  was  a  small 
variation  in  flie  inductance  gradient,  when  the  thickness  of  the  plate  was  varied,  and  a  plate  with  no 
thickness  could  be  assumed.  Althou^  these  first  two  assumptions  make  the  problem  easier  to  solve, 
they  may  introduce  an  infinite  current  density  at  some  of  the  edges.  Stickily  speaking,  this  behavior 
invalidates  the  third  assumption,  because  a  polynomial  caimot  be  infinite  at  the  edges.  Despite  this 
complication,  it  is  possible  to  get  meaningful  results  with  some  testing  and  precautions. 
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Altiiough  the  integrals  in  equation  (14)  are  in  principle  analytic,  their  explicit  expressions  quickly 
become  long  and  tedious  even  for  small  i  and  j.  In  the  end,  however,  it  is  their  values  that  are 
required  to  perform  the  calculations.  A  procedure  to  find  their  values  or  expressions  starts  by 
identifying  the  general  form  of  the  integrals  that  appear  in  all  of  diese  equations,  which  is: 


Fix,y,z)k,i 


fdx-fjy’ 

-b 


{{x-x')^  +  iy-y')^  + 


(A-1) 


Equation  (14),  for  example,  may  be  rewritten  as 

IJix,y,z)^i  =  a'^b^jyFix,y,z)ij.i  -  a^.fy(;  +  2)F(x,y,z),.^.^j 

-  b^yjF(x,y,z)i^2,j-i  yC/+2)F(x,y,z)j^2.;>i 

-  a^b^F(xy,z)ij  +  a2(;  +  l)F(x,y,z).^.^2 

+  b^jFix,y,z)i^2j  ~  0‘+2)F(x,y,z).^2j+2'  (^-2) 


Unfortunately  equation  (A-1)  is  not  the  final  integral  that  must  be  evaluated,  because  the  first 
step  taken  to  perform  the  integration  is  to  introduce  a  change  in  variables: 


F(x,y,z)i^i 


x-a  y-b 

fdufdv 

x+a  y-^b 


(x-u)^  (y-vy 


(A-3) 


where  u  =x  -  x',du  =  -dx ',  v  -y  -y',  and  dv  =  -dy Applying  the  binomial  theorem  to  the 
numerator. 


(x-  M)*(y- v)^ 


E  E  Ox*-"y'-”(-l)'‘*'"M“v'" 

„=o  m=o  V 


(A-4) 
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and  defining  a  new  set  of  integrals. 


As  an  example,  if  i  -  1  and/  =  1,  then  the  numerator  in  equation  (A-3)  becomes  (x  -  u)  (y  -  v)  = 
xy  -  yu  -  XV  +  uv  and  the  integral  becomes 

F =  xyJ{u,v\^  -yJ -xJ ^  +  J (u,v\  i  (A-7) 

in  terms  of  the  fundamental  integrals.  It  is  implied  that  these  fundamental  integrals  are  evaluated 
at  the  limits  of  integration.  In  more  explicit  terms, 

F(x,y,z\^  =  xyiJ(x-a,y-b\^-J(x-a,y*b\^-J(x+a,y-b\^*J(x+a,y*b\^) 
-y(J(x-a,y- b^^  -J{x-a,y* b\^  -J{x+a,y- b\^  +  / (x +a, y + b\  ^) 

-  xiJ  (pc- a,y-b\^-  J  (x-  a,y+b\^- J(x+a,y- b\^  +  J  (pc+a,y+b\^) 

-  (J (x- a,y- \b\  j  -  / (x- a, y  +  b\^  -J(x+a,y- b\  ^  *J(x+a,y  +  b\^).  (A-8) 

The  value  of  the  integral,  equation  (A-5),  for  any  n  and  m  ciui  be  found  by  using  a  set  of 
recursion  relations,  which  is  a  relationship  between  a  higher  order  integral  in  terms  of  lower  order 
integrals.  These  recursion  relations  start  with  the  four  lowest  order  integrals: 
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J(u,v)oo  =  -arctan  — 

Z  zr 

(A-9) 

•^(«,v)o,  =  -ln(«+r), 

(A-10) 

o 

II 

1 

5* 

+ 

(A-11) 

J  (w»v)j  j  r, 

(A-12) 

where  r  =  +  v^+z^  hi  these  e(]UEtions  snd  the  ones  following.  The  recursion  rel&tions  are  for 

n  >  0  and  m  >  2, 

(n+m-l)/(«,v)„^  =  (w-l)M”*‘V(«)^.2-nv'”-'C/(v)„-(m-l)zV(ii,v)^„.2  (A-13) 


and  for  n  >  2  and  m  >  0, 


where  U  (vX  and  V  (u)j  are  the  integrals. 


x-a  I 

Uiv).  =  j  du—  and  V{u).  = 

JC+fl 


/ 


(A-15) 


Each  of  these  integrals  have  a  recursion  relation: 

iU  (y\  =  r  -  O'  - 1)  (v^  +  z^)  U (vX._2 


and 

jV  (u).  =  v-'-^  r-(J-l)  («^  +  z^)  V  (u)j_2. 


(A-16) 

(A-17) 
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where  the  starting  integral  for  U  (v)i  are  U  (v)o  =  In  (m  +  r)  and  U  (v)i  =  r,  and  the  starting  integrals 
for  V  (M)j  are  V  {u\  =  In  (v  +  r)  and  V  (u)i  =  r.  These  recursion  relations  can  either  be  used  to  write 
an  explicit  expression  for  equation  (A-1),  or  they  can  be  used  to  write  a  computer  program,  like  the 
one  listed  in  the  next  section,  to  calculate  its  value.  This  program  does  not  calculate  the  eddy 
currents  in  the  plate,  but  it  does  demonstrate  how  these  equations  are  used. 
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#include  <stdio.h> 

#mclude  <math.b> 

/*  NMAX  is  the  absolute  maximum  value  for  nmax  */ 

#defineNMAX20 

/*  MMAX  is  die  absolute  maximum  value  for  mmax  ♦/ 

#defineMMAX20 

struct  plate  { 

double  a;  /*  The  half  width  of  the  plate  ♦/ 

double  b;  /*  The  half  height  of  the  plate  */ 

double  X,  y,  z;  /*  The  field  point  coordinates  */ 

}p; 

void  main  ()  { 

double  fityz  ( int  k,  int  el,  struct  plate  *p ); 
int  k,  el; 

/*  Ask  for  and  receive  the  half  width  */ 
printf  ("Half  width  -  "); 
scanf  ("%le",  &(p.a) ); 

/♦  Ask  for  and  receive  the  half  height  */ 

printfCHalfheigth-"); 

scanf  ("%le",  &(p.b) ); 

/*  Ask  for  and  receive  the  field  point  *! 

printf  ("Field  coordinate  (x,y,z)  -  "); 

scanf  ("%le  %le  %le",  &(p.x),  &(p.y),  &(p.z) ); 

/♦  Ask  for  and  receive  k  */ 

printf("k-"); 

scanf  ("%d",  &k ); 

/*  Test  if  k  is  too  big  */ 
if  ( k  >  NMAX  )  k  -  NMAX; 

/*  Ask  for  and  receive  1  */ 
printf  (•!-"); 
scanf  ("%d",  &el ); 

/*  Test  if  el  is  too  big  */ 
if  ( el  >  MMAX )  el  -  MMAX; 
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printf  ("fxyz  -  %e\n",  fxyz  ( k,  el,  ) ); 


}  /*  End  of  the  Main  program  */ 


/♦  This  procedure  calculates  the  integral,  equation  A1  */ 

double  fxyz  ( int  k,  int  el,  struct  plate  ♦pd ) 

{ 

void  juv  (  double  t[][MMAX],  int  n,  int  m,  double  u,  double  v,  double  z ); 
void  poly_pow  ( double  *c,  double  x,  int  n); 

double  sum; 
double  xu[NMAX]; 
double  yv[MMAX]; 
double  capj[NMAX][MMAX]; 
double  temp[NMAX][MMAX]; 

int  n,  m; 

/*  Evaluate  juv  at  the  limits  of  integration  */ 

juv  ( temp,  k,  el,  pd->x  -  pd->a,  pd->y  -  pd->b,  pd->z ); 

for  (n  -  0;  n  <-  k;  n++) 

for  (m  -  0;  m  <-  el;  m++)  capj[n][m]  -  temp[n][m]; 

juv  ( temp,  k,  el,  pd->x  +  pd->a,  pd->y  +  pd->b,  pd->z ); 
for  (n  -  0;  n  <-  k;  n-H-) 

for  (m  -  0;  m  <-  el;  m++)  capj[n][m]  +-  temp[n][m]; 

juv  ( temp,  k,  el,  pd->x  -  pd->a,  pd->y  +  pd->b,  pd->z ); 
for  (n  -  0;  n  <-  k;  n++) 

for  (m  -  0;  m  <-  el;  m++)  capj[n][m]  —  temp[n][m]; 

juv  ( temp,  k,  el,  pd->x  +  pd->a,  pd->y  -  pd->b,  pd->z ); 
for  (n  -  0;  n  <-  k;  n-H-) 

for  (m  -  0;  m  <-  el;  m-H-)  capj[n][m]  —  temp[n][m]; 

/*  capj[n][m]  matrix  now  has  the  values  for  all  the  fundamental  */ 

/*  integrals,  equation  A5,  needed  to  find  the  values  for  the  */ 

/*  integral,  equation  A3. 

♦/ 

sum  -  0.0; 

poly_pow  (  XU,  pd->x,  k);  /*  (x-uj-^k  in  equation  A3*/ 

poly_pow  (  yv,  pd->y,  el);  /♦  (y-v)^l  in  equation  A3*/ 
for  (  n  -  0;  n  <-  k;  n-H-) 
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for  (  m  -  0;  m  <-  el;  m++ )  sum  +-  xu[n]*yv[m]*capj[n][m]; 
return  sum; 

> 

/*  Performs  the  integral,  equation  A5  */ 

/*  The  results  are  stored  in  a  temporary  matrix  so  that  they  may  be  */ 

/*  added  or  subtracted  according  to  the  limit  of  integration.  */ 
void  juv  ( double  t[][MMAX],  int  k,  int  el,  double  u,  double  v,  double  z )  { 

double  capu[NMAX]; 
double  capv[MMAX]; 

double  vz,  uz,  zsq; 
double  r, 
double  num,  den; 
double  up,  vp; 
inti,j,  n,m; 

zsq«  z  *  z; 
vz  -  v*v  +  zsq; 
uz  -  u*u  +  zsq; 
r  -  sqrt  ( vz  +  u*u  ); 

/*  The  starting  values  for  the  recursion  relation,  equation  A16  */ 
if  ( u  >  0.0  )  capu[0]  -  log(u+r); 

else  capu[0]  -  log(vz/(r-u) );  /*  An  identity  used  when  u  is  negative  */ 
capu[l]  -  r, 

/*  The  recursion  relation  itself,  equation  A16  */ 
up  -  u;  /*  up  is  u  raised  to  the  i-1  power  */ 
for  (i  “  2;  i  <-  k;  i-H-)  { 

capu[i]  -  (up*r  -  (i-l)*vz*capu[i-2])/i; 
up  *-  u; 

> 

/*  The  starting  values  for  the  recursion  relation,  equation  A17  */ 
if  ( V  >  0.0  )  capv[0]  -  log(v+r); 

else  capv[0]  -  log(  uz/(r-v) );  /*An  identity  used  when  v  is  negative  ♦/ 
capv[l]  -  r, 

/*  The  recursion  relation  itself,  equation  A17  */ 
vp-v;  /*  vp  is  V  raised  to  the  j-1  power  */ 

for(j-2;j<-el;j++)  { 

capv[j]  -  (vp*r  -  (j-l)*uz*capv[j-2])/j; 
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vp  *-  v; 

> 

/*  equation  A9  */ 
num  -  u*v; 
den  -  z  ♦  r; 

/*  Test  if  z  is  close  to  zero  ♦/ 

if  ( fabs  ( z )  <  l.Oe-32 )  t[0][0]  -  0.0;  /♦  The  limiting  value  */ 
else  t[0]  [0]  -  atan2  (num,  den)  /  z; 

/*  equation  AlO  */ 

if  ( u  >  0.0 )  t[0][l]  -  -log(u+r); 

else  t[0][l]  -  log((r-u)/vz)*y*  An  identity  used  whra  u  is  negative  */ 

/*  equation  All  */ 

if  ( V  >  0.0 )  t[l][0]  -  -log(v+r); 

else  t[l][0]  -  log((r-v)/uz)y*  An  idaitity  used  when  v  is  negative  */ 

/♦  equation  A12  */ 
t[l][l]  -  -r. 

/*  Use  equation  A13  to  fill  out  the  top  row  of  the  matrix  ♦/ 
n-0; 

for  ( m  -  2;  m  <-  el;  m++ ) 

t[n][m]  -  u*capv[m-2]  -  zsq*t[n][m-2]; 

/*  Use  equation  A13  to  fill  out  the  second  row  from  the  top  */ 
n-  1; 
vp-v; 
up  -  u*u; 

for  ( m  -  2;  m  <-  el;  m++ )  { 

t[n][m]  -  ((m-l)*up*capv[m-2]-vp*capu[n]-(m-l)*zsq*t[n][m-2])/m; 
vp  *-  v; 

> 

/*  Use  equation  A14  to  fill  out  each  column  of  the  matrix  ♦/ 
vp-v; 

for  (  m  -  0;  m  <-  el;  m++  )  { 
up-u; 

for  ( n  -  2;  n  <-  k;  n++ )  { 

t[n][m]  -  (n-l)*vp*capu[n-2]-m*up*capv[m]-(n-l)*zsq*t[n-2][m]; 
t[n][m] /- n+m-1; 
up  *-  u; 

> 

vp  *-  v; 

> 

)/*Endofjuv*/ 
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/*  (x-u)  raised  to  the  n-th  power  for  a  given  x  value  */ 

void  poly_j)ow  ( double  *p,  double  x,  int  n ) 

{ 

/♦  The  index  of  p[i]  is  the  power  of  u  or  v;  p[i]  means  u^i  or  v^i  */ 

/*  The  value  of  p[i]  is  the  coefficient  for  the  term  */ 

/*  Example:  p[0]  -  1.0,  p[l]  -  -2.0,  and  p[2]  -  1.0,  would  represent  */ 
/*  the  polynomial  1.0  -  2.0*u  +  1.0*u*u  or  (l-u)^2  for  x  -  1  */ 

inti; 

if(n  — 0){ 

p[0]-1.0; 

return; 

> 

/*  -1  raised  to  the  n  th  power  */ 
if  (0x0001  &  n)p[n]--1.0; 
else  p[n]  - 1.0; 

for  ( i  -  n-1;  i  >-  0;  i--)  p[i]  -  -x  *  p[i+l]  *  (i+1)  /  (n-i); 

}  /*  End  of  poly_pow  ♦/ 
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